dists <- readRDS("temp/dists_long_new_from_gva.rds") |> 
  mutate(year = year(as.Date(date)),
         id = id4, ## id (and distance, in next line) of nearest 4-victim shooting
         dist = dist4) |> 
  filter(!is.na(id),
         id %in% (readRDS("temp/geocoded_shootings_from_gva.rds") |> 
                    filter(home + dv == 0))$id, ## drop if in the home or domestic violence,
         dist4 <= 10
  ) |> 
  group_by(GEOID, year) |> 
  tally()


bg_share_dem_16 <- function(s){
  library(tigris)
  library(rgdal)
  library(sf)
  library(sp)
  library(rgeos)
  library(SearchTrees)
  library(raster)
  library(data.table)
  library(tidyverse)
  
  state <- substring(list.files(s)[1], 1, 2)
  
  if(!file.exists(paste0("temp/", state, "_bg_prec_16.rds"))){
    scode <- fips_codes |>
      rename(s = state) |> 
      filter(tolower(s) == tolower(state)) |> 
      select(state_code) |> 
      distinct() |> 
      pull()
    
    
    bgs_s <- filter(dists, substring(GEOID, 1, 2) == scode, year == 2016)
    
    if(nrow(bgs_s) > 0){
      bg_shapes <- block_groups(state = state, class = "sf", year = 2019)
      bg_shapes <- subset(bg_shapes, GEOID %in% bgs_s$GEOID)
      
      centroids <- st_centroid(bg_shapes)
      
      ps <- read_sf(dsn = s, layer = substring(list.files(s)[1], 1, nchar(list.files(s)[1]) - 4))
      colnames(ps)[1:(ncol(ps) - 1)] <- toupper(colnames(ps))[1:(ncol(ps) - 1)]
      
      ps <- st_transform(ps, crs = st_crs(centroids)) |> 
        mutate(share_dem = G16PREDCLI / (G16PRERTRU + G16PREDCLI)) |> 
        select(share_dem)
      
      tryCatch(
        {
          ps <- st_make_valid(ps)
        }, error = function(e) {
          Sys.sleep(1)
          cat("Error:", e$message, "\n")
        })
        
      ps <- ps[st_is_valid(ps), ]
      
      bg_shapes <- st_join(centroids, ps, join = st_within) |> 
        select(GEOID, share_dem) |> 
        st_drop_geometry() |> 
        mutate(year = 2016)
    }else{
      bg_shapes <- data.frame(GEOID = NA, share_dem = NA, year = NA)
    }
    
    saveRDS(bg_shapes, paste0("temp/", state, "_bg_prec_16.rds"))
  }else{
    bg_shapes <- readRDS(paste0("temp/", state, "_bg_prec_16.rds"))
  }
  
  return(bg_shapes)
}

bg_share_dem_20 <- function(s){
  library(tigris)
  library(rgdal)
  library(sf)
  library(sp)
  library(rgeos)
  library(SearchTrees)
  library(raster)
  library(data.table)
  library(tidyverse)
  
  state <- substring(list.files(s)[1], 1, 2)
  
  if(!file.exists(paste0("temp/", state, "_bg_prec_20.rds"))){
    scode <- fips_codes |>
      rename(s = state) |> 
      filter(tolower(s) == tolower(state)) |> 
      select(state_code) |> 
      distinct() |> 
      pull()
    
    
    bgs_s <- filter(dists, substring(GEOID, 1, 2) == scode, year == 2020)
    
    if(nrow(bgs_s) > 0){
      bg_shapes <- block_groups(state = state, class = "sf", year = 2019)
      bg_shapes <- subset(bg_shapes, GEOID %in% bgs_s$GEOID)
      
      centroids <- st_centroid(bg_shapes)
      
      ps <- read_sf(dsn = s, layer = substring(list.files(s)[1], 1, nchar(list.files(s)[1]) - 4))
      colnames(ps)[1:(ncol(ps) - 1)] <- toupper(colnames(ps))[1:(ncol(ps) - 1)]
      
      ps <- st_transform(ps, crs = st_crs(centroids)) |> 
        mutate(share_dem = G20PREDBID / (G20PRERTRU + G20PREDBID)) |> 
        select(share_dem)
      
      tryCatch(
        {
          ps <- st_make_valid(ps)
        }, error = function(e) {
          Sys.sleep(1)
          cat("Error:", e$message, "\n")
        })
      
      ps <- ps[st_is_valid(ps), ]
      
      bg_shapes <- st_join(centroids, ps, join = st_within) |> 
        select(GEOID, share_dem) |> 
        st_drop_geometry() |> 
        mutate(year = 2020)
    }else{
      bg_shapes <- data.frame(GEOID = NA, share_dem = NA, year = NA)
    }
    
    saveRDS(bg_shapes, paste0("temp/", state, "_bg_prec_20.rds"))
  }else{
    bg_shapes <- readRDS(paste0("temp/", state, "_bg_prec_20.rds"))
  }
  
  return(bg_shapes)
}


cl <- makeCluster(8)  
registerDoParallel(cl)

clusterExport(cl, list("bg_share_dem_16", "dists", "bg_share_dem_20", "bg_share_dem_16_PRE"))

runs <- list.dirs("../raw_data/vest/vest_2016/2016")[2:length(list.dirs("../raw_data/vest/vest_2016/2016"))]

out <- rbindlist(parLapply(cl, runs,
            fun = bg_share_dem_16))

runs <- list.dirs("../raw_data/vest/vest_2020")[2:length(list.dirs("../raw_data/vest/vest_2020"))][c(1:32, 34:52)]

out <- rbindlist(parLapply(cl, runs,
                           fun = bg_share_dem_20))